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ABSTRACT 

We model the comet ary structure around Mira as the interaction of an AGB wind from Mir a A, and 
a streaming environment. Our simulations introduce the following new element: we assume that after 
200 kyr of evolution in a dense environment Mira entered the Local Bubble (low density coronal gas). 
As Mira enters the bubble, the head of the comet expands quite rapidly, while the tail remains well 
collimated for a > 100 kyr timescale. The result is a broad- head/narrow-tail structure that resembles 
the observed morphology of Mira's comet. The simulations were carried out with our new adaptive 
grid code Walicxe, which is described in detail. 

Subject headings: circumstellar matter — ISM: jets and outflows — hydrodynamics — methods: 
numerical — stars: AGB and post AGB — stars: individual (Mira) 
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1. INTRODUCTION 

Mira is a well studied binary system. Mira A is a proto- 
typical, thermally pulsating luminous AGB star (which 
in fact gives name to the so-called Mira class of vari- 
ables), while Mira B is a much less luminous star, be- 
lieved to be a white dwarf or a main sequence star. Re- 
cently, the discovery of a cometary structure of ~ 2° 
in the sky w ith its head centered on Mira, by |Martin 
et al. (2007), with the GALEX satellite, has generated 
a renewed interest in this system. In order to compare 
the observations with our models we inclu de in Figure [T 
a sche matic diagram of the UV image of Mar tin et al. 
([2007). 



The broad-band UV filters (FWHM AA = 256, 730 A, 
centered at A c = 1516, 2267 A) in GALEX make difficult 
to discriminate the phy sical mechanism resp onsible for 
the observed emission. Martin et al.| ( |2007| ) suggested 
that most of the observed emission in the tail of Mira 
could be due to fluorescent H2 lines, while in regions 
closer to the bow-shock it is likely that other species, 
such as C IV could contribute to the emission. 

The distance to Mira, estimated from Hipparcos is 
107 pc ( Knapp et al.|2003| . The system Mira moves at a 



high speed (~ 125 km s ) with respect to its sur round 



ing medium , as deduced from i ts proper mo tion (Turon 
et al-|1993b and radial velocity flEvans||1967| ). The phys 



ical size of the cometary tail, assuming the Hipparcos 
distance is ~ 4pc. 

New observations and/or reexamination of previous 
observations have been made in several wav elenghts, in - 
cluding H I 21cm (Ma tthews et al.|2008L IR flUeta|2008| ), 
and optical ( jMeapurn et al. 2009)' Recen t theoretical 
work h as been done by Wareing et al. (2007); Raga et al 



(|2008j); |Raga fc Cant6| ( |2008p .~rWareing et a. 



,ga et al. 
. (2007) 
the in- 



performed a 3D hydro dynamical simulation of 
teraction between an isotropic AGB wind and the ISM, 
modeled as a plane-parallel wind (considering a reference 
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system in which Mira is at rest). Raga et al. (2008) in- 
cluded an inhomogeneous AGB wind, resulting m a more 
complicated shock structure for the cometary head (pos- 
sibly more consistent with the observations). These two 
models, however, have one flagrant flaw: if the models 
are adjusted to produce the correct size for the cometary 
head, the tail that they produce is significantly broader 
than the one of Mira's comet. 

The reason for this discrepancy with the observations 
lies in the fact that the ISM density in those models 
is estimated by assuming ram pressure balance between 
Mira's wind and the ISM at the tip of the head of the 
cometary structure. Therefore, the distance between the 
front of the bow-shock and Mira is consistent with the 
observations. However, the tail produced by the model 
has a width comparable to the head, while Mira's comet 
s hows a considerably n arrower tail (see Fig [I]). 



Wareing et al. (2007) speculated that Mira could have 
recently crossed into the Local Bubble (a tenuous region 
of coronal gas). Most of the tail could then have formed 
outside the Local Bubble, in a higher density ISM, result- 
ing in the observed, narrow tail structure. In the present 
paper, we explore this possibility through a set of axisym- 
metric numerical simulations of a travelling star+wind 
system which emerges into a hot, coronal gas bubble. 

The paper is organized as follows: in section 2 we de- 
scribe the model, including a review of the analytical 
solutio n of the wind / st r eamin g environment by Wilkin 
(1996) and Canto et al. (1996). The setup of the numer- 
ical simulations is presented in section 3, followed by the 
results in section 4, and a summary in section 5. As this 
paper is the first publication with our new hydrodynam- 
ical code Walicxe, we include an appendix describing 
the code in detail, as well as some tests. 

2. THE MODEL 

Following [Wareing et al] fl2007| ) and |Raga et al.| ( |2008 ) 
we model Mira's cometary tail as the interaction of a 
slow, dense, wind that arises from the AGB star Mira 
A, with a fast, lower-density, plane-parallel wind that 
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Fig. 1. — Sketch of the UV image of the tail of Mira observed with GALEX, fig. la of Mar tin et al.| (2007). The intensity is represented 
by the different shades of gray (darker corresponds to a higher intensity). The orientation (!N,hJ) and proper motion (PM) are indicated by 
the arrows below the head of the comet. 



iow-shock 




Streaming ISM 
A 



Fig. 2. — Schematic of the two- wind interaction model. The 
interaction of a moving spherical wind source (Mira A) and its 
surrounding ISM can be modeled, in a reference frame with Mira 
A at rest, as a stationary wind source accompanied by a streaming 
ISM. In this figure, we also illustrate the spherical radius R used in 
the analytical model, the polar angle 0, and the standoff distance 
Rs- 

results from the motion of Mira through the ISM (the 
simulations are performed in the reference frame in which 
Mira is at rest), as shown in the schematic diagram of 
Figure |2| 

The structure that results from the interaction between 
a streaming, plane-parallel flow and a spherical wind has 
a cometary shape that has been described analytically, 
considering conservation of linear a nd angular momen- 
tum in a "thin-shell" approximation ( |Wilkin|1996 Canto 



et al. 1996). With this formalism, the shape of bow- shock 



is found to have the form: 



R(0) = R s esc (9^3(1 -0cot0), 



(1) 



where R is the spherical radius (measured from the cen- 
ter of the isotropic wind source), 6 is the polar angle 
(measured from the symmetry axis, aligned with the di- 
rection towards the impinging, plane-parallel flow), and 
R s is the stagnation radius (or standoff distance, see Fig- 
ure]^. The stagnation radius (point of closest approach 
to the spherical wind source) is given by the ram-pressure 
balance between the two winds: 



(2) 



where M w is the mass loss rate of Mira, v w its terminal 
wind velocity, p a the ambient (ISM) mass density, and v a 
the velocity of the star moving through the surrounding 
environment. 




Fig. 3. — Analytical solution of t he two-wind "thin-shell" model 
(Wilkin) [19961 [Canto et al] [1996b for three different stagnation 
radii: R s = 0.1, U.5, l.U (solicT, clotted and dashed lines, respec- 
tively). The position of the spherical wind source is at (0,0), 
marked by the star. The plane-parallel wind impinges form the 
right as in Fig. [2] 

Figure [3] shows the locus of the interaction region (lo- 
cation ofl the thin-shell, and shape of the bow-shock 
structure) for different stagnation radii. Notice that for 
smaller R s the cometary tail becomes narrower. 

3. THE NUMERICAL SETUP 

We performed two hydrodynamical simulations of 
Mira's tail with the new code Walicxe. The code solves 
the gas-dynamic equations in a block based adaptive 
mesh that is designed to run in parallel on machines with 
distributed memory (i.e. clusters). 

Walicxe can be run with a Cartesian two-dimensional 
mesh, or an axisymmetric (cylindrical) grid. For this par- 
ticular problem we have used the axisymmetric version 
of the code with the hybrid HLL-HLLC Riemann solver 
(see the Appendix for details). 

Along with the gas-dynamic equations, a rate equa- 
tion for hydrogen is solved. The resulting H ionization 
fraction is used to compute the radiative energy losses, 
which are obtained with the parametrized cooling func- 
tion (that depends on the temperatur e, density and hy- 
drogen ionization fraction) described in Raga & Reipurth 
( [2004] ). 
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The simulations have 5 root blocks of 16 x 16 cells 
each, aligned along the axial direction (z— axis). We al- 
low 7 levels of refinement, yielding an equivalent reso- 
lution of 1024 x 5120 (radial x axial) cells at the high- 
est grid resolution. The computational domain extends 
[(0, 0.25) x (0, 1.25)] x 10 19 cm in the radial and axial di- 
rections, respectively. Thus, the resolution at the highest 
level of refinement is « 2.4 x 10 15 cm. We impose a reflec- 
tive boundary condition along the symmetry axis (r = 0) . 
Open boundary conditions are used at r = 0.25 x 10 19 cm 
and z = 0. 

In both models the spherical wind has the following 
properties. It is injected (at every timestep) in a spher- 
ical region of 3 x 10 16 cm radius, centered at r = 0, 
z = 1.05 x 10 19 cm. The density inside the injection 
region follows an oc R~ 2 law (R is the spherical radius 
measured from the position of the wind source) , scaled so 
that the mass- loss rate is M w = 3 x 10 -7 M yr _1 . The 
wind is neutral, has a temperature of T w = 10 2 K, and 
an outward, radially directed velocity v w = 5 km s _1 . 
The values of the physical parameters of Mira's wind 
were taken fr om t he observatio nal studies of CO lines by 



Y oung fll995| ) and |Ryde et al.| ( |2000D 

At the initial time (t = 0), both models are identical, 
with a plane-parallel wind that fills the computational 
domain (except inside the sphere in which the isotropic 
wind is imposed), which is replenished by an inflow con- 
dition at z = 1.25 x 10 19 cm. The ISM is homogeneous, 
with a hydrogen number density nn = 1 cm -3 (consis- 
tent with the ISM in the galactic plane, j ust outside the 
Local Bubble in the vicinity of Mira, see Ver gely et al. 
200l[ ), and a velocity v a = 125 km s _1 (in the — z direc- 
tion). The wind and ISM are neutral, except for a seed 
electron density assumed to arise from singly ionized C 
(i.e. the minimum ionization fraction is 10 -4 ). The ISM 
has a temperature of T a = 10 3 K. 

In model Ml, we evolve these initial/boundary condi- 
tions up to an integration time of 350 kyr. Model Ml 
corresponds to the evolution of Mira outside of the Local 
Bubble, in a uniform medium dense enough to produce 
the narrow tail as observed. 

In model M2, we simulate Mira entering the Local Bub- 
ble, a region much hotter and tenuous than the aver- 
age ISM. To achieve this, after 200 kyr of evolution in 
the same medium as model Ml, we raise the tempera- 
ture of the ISM that enters the computational domain to 
T a = 10 6 K, lower its den sity to nn = 5 x 10~ 2 cm 



(from the observations by Lallement et al. 2003), and 
change its ionization fraction to 1 (fully ionized nydro- 
gen). We then let this configuration evolve up to a 
t = 350 kyr integration time. 

The model of Mira entering the Local Bubble is well 
justified by the current Galactic Position of Mira (l, b) = 
(168,-58) flPerryman et al.l|1997| [Wareing et al.||2QQ7| ). 
This position puts Mira inside, but close to the edge of 
th e Local Bubble, as can be seen from the maps presented 



by Lallement et al. (2003). Given the high velocity of 
Mira, in the ~ 300 kyr that it takes to form a tail of 
4 pc, it would have traveled approximately 38 pc. This 
distance is enough to place the starting point of the tail 
formation outside the Local Bubble. 



4. RESULTS 



We let the models to evolve up to an integration time of 
350 kyr, which is of the order of the time needed to form 
a cometary tail of ~ 4 pc in length (the size of the struc- 
ture observed with GALEX). Figure [4] shows the density 
(top half of the panels) and temperature distributions 
(bottom half) of Ml and M2 at the final evolutionary 
stage. 

From the top panel in Figure [4] we can see that after 
350 kyr a narrow tail, of axial and radial dimensions 
comparable with the observations of Mira, is formed. 
However, the stand-off distance (in agreement with the 
analytical prediction in eq. [2| is R s ~ 4.9 x 10 16 cm, 
which is much closer to the source compared to the 
R s ~ 3.1 x 10 1T cm that can be deduced from the obser- 
vations of Mira's comet. The analytical solution for the 
parameters of model Ml, plotted as a dotted line, traces 
well the region that separates the dense and cold mate- 
rial which arises from Mira's wind, and the hotter, lower 
density shocked ISM. This can be seen more clearly in 
the following section, where we separate the contribution 
of the wind and show column density maps. 

The bottom panel of Figure [4] shows the tail that 
would be produced if Mira entered the Local Bubble 
some 145 kyr ago (i. e., model M2). The analytical so- 
lution for the parameters used for the Local Bubble is 
included as a dashed line. The corresponding standoff 
distance for model M2 is at R s = 2.2 x 10 17 cm to the 
right of Mira, a distance that is close to what is observed 
(3.1 x 10 17 cm). 

This standoff distance, and the corresponding width of 
the tail (dashed line) are much larger than observations. 
However, the change from one stationary configuration 
to the other is not instantaneous. As the cometary flow 
expands after entering the Local Bubble the head grows 
rapidly. Approximately 20 kyr after the entrance of Mira 
into the Local Bubble we obtain a structure that re- 
sembles remarkably well the broad- he ad /narrow- tail of 
Mira's comet, such tail remains well collimated for at 
least 100 kyr. It is only after 150 kyr of entering the Lo- 
cal Bubble that the tail becomes somewhat wider than 
the UV observations (see time sequence in the following 
section). 

A clear difference in the wind/environment interface 
can be seen in models Ml and M2. In model M2, the 
outer boundary of the region filled by the stellar wind 
shows a more complex structure, which fragments the 
tail downstream from the wind source. Such structures 
are absent in Model Ml (see Figure [4]). This difference 
is a result of the fact that while the wind/environment 
boundary has a highly supersonic velocity shear in model 
Ml, this shear is barely transonic in model M2 (as Mira's 
motion is basically sonic with respect to the hot, Local 
Bubble environment). It is a well known experimental 
result that the spreading rate of sonic mixing layers is 
much larg er than the one of h igh Mach number mixing 
layers (see |Canto fc Ra ga 1991), this effect is also present 
in our models, therefore a more turbulent tail is seen in 
model M2. 

4.1. Column density maps 

One problem of interpreting the |Martin et al.l ( |2007[ ) 
observations of Mira is that the broad band filters m 
GALEX do not give much information about the physical 
mechanisms that produce the observed emission. They 
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(a) M1, t = 350.0 kyr 
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Fig. 4. — Density (top half of each panel), and temperature (bottom half) stratifications at at = 350 kyr integration time, for the two 
models. The numerical values are shown in the logarithmic gray-scale (color-scale in the on-line version) bars at the right of each panel. 
We have included in both panels the analytical solution for the bow-shock for the two stagnation radii used (see text for details): the dotted 
line is for R s = 4.9 x 10 16 cm and the dashed line for R s = 2.2 x 10 17 cm. The temperature panel has been reflected on the symmetry axis 
for visual purposes. 



suggested that the bulk of the emission arises from H2 
molecules, which are likely to be present in the cold wind 
of Mir a, but not in the surrounding environment. 

In order to make a better comparison with the obser- 
vations we have added a passive scalar in the simulations 
that allows us to identify the material from Mira's wind 
and the ISM material. We then computed column den- 
sity maps of the material from Mira's wind (which would 
be molecular H2) and present them in the form of a time 
sequence in Figure [5] 

In Figure [5] we see that the tail forms and reaches the 
observed length (of Mira's comet) after approximately 
~ 300 kyr. The ISM that enters from the right changes 
(in model M2) to the Local Bubble parameters at t = 
200 kyr, but it takes some ~ 5 kyr to reach the two- wind 
interaction region. From then on, models Ml and M2 
start to differ. 

From the column density maps it is clear that the stel- 



lar wind material in Ml is well confined to a region that 
resembles the tail of Mira for the duration of the simu- 
lation, but (as we have discussed above) this requires an 
ambient density that puts the axial standoff distance too 
close to Mira. 

For model M2, as Mira crosses into the Local Bubble its 
wind expands towards a new steady state configuration 
(which is not reached in our simulations). We found that 
at after only 20kyr the stellar wind has expanded enough 
to form a dense structure extending ~2x 10 17 cm up- 
stream from the stellar wind source. This size is in agree- 
ment with the axial extend of the head of Mira's comet. 
We also see that the tail in model M2 expands at a much 
slower rate compared to the head, therefore reproduc- 
ing the "broad head/narrow tail" morphology of Mira's 
comet for a timescale on the order of 100 kyr. After 
150 kyr inside the Local Bubble the tail of Mira has ex- 
panded laterally beyond what is observed in the GALEX 
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Fig. 5. — Time sequence of H2 column density (i. e., of the material from Mira's wind), for models Ml and M2. The value of the column 
density can be read from the logarithmic gray-scale bar (color-scale in the on-line version) at the top of the plots. The time of evolution is 
indicated in the legend at the upper left corner of each panel. For the first three panels (a-cf) the evolution is identical for both models. We 
have included a blow-up of the region close to the wind source for an integration time of 300 kyr (/and g). The physical size is indicated 
by the arrows in panels (f) and (g). 



UV maps (see Martin et al. 2007, and diagram in Fig. 

0. 

To simulate the UV emission that would arise from 
the H2 molecules requires a chemical/radiation transfer 
calculation that is beyond the scope of this paper. We 
only note that if we use the standard conversion factor 
Ay « 10 -21 7Vh (due to dust absorption), we find that 
the tail in our simulations has Ay ~ 10. Therefore, 
molecules present in the wind from Mira A would not 
be photodissociated in a substantial way, and the tail 
should indeed have a large H2 fraction. 

Two interesting issues arise from this work, one is the 
age of Mira's tail itself, the other is how long it has been 



inside the Local Bubble (assuming that this is the case). 

Regarding the age of the tail, in the original discov- 
ery of |Martin et al. ( |2007 ), the authors estimate that 
the tail was formed in apprc 



appro ximately 30 kyr, while the 



simulations of Wareing et al. ( |2007[ ) suggested an age of 
450 kyr. This large discrepancy is related to the rate 
at which Mir a's wind dece l erates to merge dynamically 
with the ISM. |Martin et alT ( 2007[ ) assumed that the wind 
material instantly decelerates, so the 30 kyr is the time 
that Mira takes to traverse th e length of the tail. In con- 
trast, the 450 kyr timescale of Wareing et al. (2007) was 
directly obtained from numerical simulations which are, 
however, sensitive to viscous effects. Such effects depend 
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mainly on two things: the density of the ISM that is 
used (a physical effect), and the amount of artificial vis- 
cosity (a numerical effect) needed to stabilize the code. 
Artificial viscosity usually overwhelms real viscosity, re- 
sulting in an enhanced deceleration, thus the numerical 
estimates correspond to lower limits of the actual age. 
Our simulations lie in between these two previous esti- 
mates, as they require a time of ~ 300 kyr to develop 
Mira's tail. This difference with the model of Waering et 
al. (2007) is due to the higher ISM density of our models, 
which results in a higher deceleration rate. 

The question of the time spent by Mira inside the Local 
Bubble is not well constrained by our models because 
the tail takes on the order of 100 kyr to expand laterally 
in an appreciable way. Our results suggest that Mira 
could have crossed into the Local Bubble between 20 and 
120 kyr ago. This result is in rough agreement with the 
130-225 kyr proposed by |Wareing et al.| ( [2007| . In this 
respect, a more detailed map of the Local Bubble and an 
extrapolation of Mira's motion might prove more useful. 

5. SUMMARY 

We have modeled the cometary structure observed 
around Mira as the interaction of the cold, dense wind 
from Mira A with a streaming ISM. For this purpose we 
have used a newly developed code, which is described in 
detail in the Ap pendix. Similar simulations have been 
done previously ( | Wareing et al.|[2007 iRaga et al. 2008), 
but they fail to reproduce the broad- head/ narrow-tail 
configuration characteristic of Mira's comet. 

The new element in our simulations is to assume that 
Mira has recently crossed (from a denser than previ- 
ously considered) medium into the Local Bubble, a re- 



gion of tenu ous coronal gas . Thi s possibility was first 
suggested by Wareing et al. (2007), and is quite feasible 
given Mira's current location. 

By starting a stellar wind in a denser region than pre- 
vious papers we are able to reproduce the width of the 
observed tail of Mira's comet. By allowing the stellar 
wind source to enter (at a later time) a less dense region 
of the ISM (i. e., the Local Bubble) we obtain an expan- 
sion of the cometary head. For the parameters of our 
model, the predicted cometary structure has a striking 
resemblance to Mira's comet after ~ 20 kyr time after 
entering the Local Bubble, and the structure last for an- 
other ~ 100 kyr. 

The age that we obtain for the formation of the tail 
is ~ 3QQ kyr. This can be compared with estima tes by 
Martin et al"] ( [20071 30 k y r ) and by | Wareing eTal] ( [20071 
450 kyr), that assume either total, or very little (respec- 
tively) deceleration of the wind as it encoun ters the ISM. 

We should note that Meaburn et al. ( 2009 ) recently dis- 
covered a bipolar jet system in the Mira system (possibly 
ejected from Mira B). This collimated outflow might have 
an important effect in the formation of Mira's cometary 
head and tail, which has not been considered in our 
present work, and should clearly be studied in the fu- 
ture. 
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APPENDIX 

A. THE WALICXE CODE 

Our previous adaptive grid codes (e. g., the YGUAZU-A, see Raga et al. 2000) were not suited for running in parallel 
computers with distributed memory (i.e. clusters). In order to achieve high resolution by today's standards we have 
developed the new code walicxe, 1 which contains the cooling/chemistry modules of YGUAZU-A, but uses a new 
adaptive mesh designed to be easily parallelizable in clusters. At this point we only have a 2D version, with a 3D 
version in progress. We describe below the numerical algorithm and the AMR scheme of the walicxe code. 



= S, 



A.l. Basic equations 
The 2D Euler system of equations can be written in Cartesian coordinates as: 

<9U <9F <9G 
dt dx dy 

where 
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(A2) 



U are the so-called conserved variables, p is the mass density, (ix, v) are the velocity components in the (x, ^-directions, 
ni, n2, n r are densities of additional species that are advected passively into the fluid, which can be used to compute 

1 The name is one of the few remaining words of the language of the charruas, a South- American tribe than inhabited what is now 
Uruguay and part of Argentina and Brazil. The meaning is "sorcery", which fits to some extent the job description of theoretical/numerical 
astrophysicists. 
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heating/cooling rates. F and G are the fluxes in the x and y directions, respectively; 

P = (m + n e ) kT 
is the thermal pressure, and the (total) energy is given by 



E ■ 



1 



P(u 2 



C V P, 



(A3) 



(A4) 



where C v is the specific heat at constant volume. S is the source vector, that contains the energy gains and losses (G 
and L, from radiative and collisional processes) and the reactions rates for the ni,ri2, ...n r species, the source vector 
can be also extended to include geometrical terms for cylindrical and spherical coordinate systems. In addition, one 
can define a vector with the so-called primitive variables 



W = 
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(A5) 



which are easily obtained from the conserved variables. 

A. 2. Discretization: finite- areas 

We will denote the indices and grid spacings of the computational cells in the x, ^/-directions by z, j and Ax, Ay, 
respectively. The time is discretized with an index n such that At n = t n +i — t n , which in principle varies at each 
time-step, however to ma ke t he notation clearer we will drop the sub-index n from At n . 

Integration of equation Al on the area of the computational cell centered at i , j, on a time interval [t + At] yields 
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the following exact expression: 
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where the variables U^- are now area averages at a time £ 

/-2/+1/2 fX+1/2 

u n = / / 

hJ Ax Ay Jy-i/2 Jx-i/2 

where x ± 1/2 = Xi ± Ax/2, y ±1/2 = yi± Ay/ 2 are the locus of the boundaries of the cell at x 



+ 1/2 



U(x,y,t n ) dxdy, 



now time (and length) averages: 
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The fluxes are 
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And the source terms are time and area averages: 
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(AlO) 
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The method we have described if often referred to as finite-vol umes , but s ince we are dealing with 2D it is more 
accurate to call it finite-areas. Of course, the expressions of eqs. (A6)-(A10) are easily extended to three dimensions 
yielding a truly finite- volumes method, where the area integrals used above become volume integrals, and the length 
integrals become surface integrals. What remains now, is to find approximations to the numerical intercell fluxes 



71+1/2 



G r 



71+1/2 



: i+1/2 j' Kjr ij+i/2i an d sources (S^ 1 ^ 2 ). This can be done by solving the Riemann problem at the cell interfaces, 
that is, assuming a true discontinuity of the hydrodynamic variables at such interfaces ( Toro|[l 999). 

A. 3. Solution method: Finite areas, second order Godunov scheme 

We use a second order Gudunov type method to advance the solution, it uses a linear reconstruction (with slope 
limiters) of the primitive variables and an approximate Riemann solver to compute the fluxes at the intercell boundaries. 
As it is customary, below we explain the method in one dimension, the expressions for the fluxes in the ^/-direction are 
analogous to those in the x-direction. The finite-volume (finite-lengths in ID) solution is: 

At / F n+l/2 _ F n+l/2 
\ r i+1/2 r i-1/2 
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To obtain a second order approximation of the intercell fluxes we proceed a s follows. First, we c alculate the timestep 
At in order to ensure that the standard Courant-Friedrichs-Lewy condition ( |Courant et al.|[l 967) is met. Then, a first 
order half-timestep is computed as: 



TT n+l/2 _ jjn ^ /™ 



At 



Sn 



(A12) 



where the fluxes are obtained solving the Riemann problem with an approximate Riemann solver (see 
using the primitive variables W^(U^). That is 

Fi — Riemann (W i} W i+ i) . 



Toro 



1999) 



(A13) 



With the values of U™ +1/2 from eq. (A12) we obtain a new set of primitives W™ +1 ^ 2 . We then use a linear recon 
struction to extrapolate them to the cell boundaries 
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(A14a) 
(A14b) 



where avg is an averaging function. The choice of this function is made in order to limit the slope used and avoid 
spurious oscillations. The code has a number of different averaging functions, the most widely known (although the 
most diffusive) is the minmod func tion. 



The left and right states in eqs. (A14) are used to approximate the fluxes used in eq. (All): 
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(A15a) 
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At this point a new At is calculated, and we iterate until the desired evolutionary stage is reached. The time step is 
the same for all blocks, independently of their level of resolution, this has a cost in computing time, and introduces 
additional numerical viscosity to coarser blocks (because they are effectively run at a smaller Courant number), but it 
makes easier and more efficient to balance the load among different processors. 

As to the Riemann solver used , the default is the HLLC algorithm, an improvement proposed by |Toro et al.| ( 1994) 
to the HLL (Harten et al. 1983) solver, which restores the information of the contact wave into the solver (the C 
stands for contact). The original HLL is also available in the code, but is rarely used because it is too diffusive. In 
addition to this solvers we have implemented a hybrid scheme that combines the HLL and HLLC fluxes to fix the 
so-called carbuncle or shock instability problem, which occur s when strong shocks align with the com putational grid. 
Our HLL-HLLC scheme is v ery s imilar to those proposed by Kim et al. (2009); Huang et al. (2010), the description 
and some tests are given in §A.6| . 

A. 4. The adaptive grid 

The hydrodynamic+rate equations described above are integrated on a binary, block-based computational grid with 
the following characteristics. 

A.4.1. Blocks (patches) 

The adaptive mesh was designed to be used parallel computers with distributed memory (i.e. clusters). It consists 
of a series of "root" blocks with a predetermined (set at compile time) number of cells, n x x n y . Each block can be 
subdivided in a binary fashion into four blocks with the same number of cells of the parent block. We will refer to 
these new blocks as "siblings" , which have twice the resolution of their parent block. The maximum number of levels 
of refinement allowed (ni evs ) is set at compile time. Thus the maximum resolution is = Z^/[2 nZeus-1 (r^)], where Li 
is the size of the computational mesh on the ith direction (i.e. i = x, y). In addition, each block has its own ghost 
cells, two in each direction for the second order scheme described below, but can be changed to implement higher order 
schemes. These cells would be communicated with MPI when adjacent blocks are in different processors. 

The result is a tree structure (of blocks or patches, not individual cells), the root blocks branch into higher resolutions 
up to ni evs . The equations are only solved in "leaf" blocks, that is, only on blocks that are not further refined. 

A. 4. 2. Refinement/ coarsening criteria 

Before advancing the solution, at every timestep we update the mesh. This is done by sweeping all the leaf blocks, 
measuring pressure and density gradients on each of them. If either the density or pressure suffer a relative gradient, 
between two adjacent cells, above a threshold (user defined, typically of ~ 0.05) the block is flagged for refinement. In 
addition, we enforce that neighboring blocks differ at most on one level of refinement, blocks in proximity to a higher 
resolution that has been marked for further refinement are flagged as well. Coarsening is also controlled by measuring 
gradients of density and pressure, but with a different tolerance (also user defined, typically of ~ 1.0), if all the relative 
gradients in the block are below this value the block is marked to be coarsen. Coarsening is only allowed when an 
entire group of four siblings have been marked and it is not impeded by the proximity criterion. 
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The actual refinement algorithm assigns the same value of the variables in the parent block to all the siblings 
(Oth— order interpolation), whereas the coarsening assigns the mean value of each variable from the four siblings to 
their parent block. 

Additional refinement criteria are easy to implement depending of the problem to be solved. One can, for instance, 
use a vorticity criteria, refine above a certain value of any variable, in certain fixed regions, or in regions labeled by a 
passive scalar advected into the flow. 

A. 4. 3. Parallelization and load balancing 

Since each of the blocks have the same number of cells, it is straightforward to send different blocks (actually a 
group of blocks) to different processors. This is achieved with the standard Message Passing Interface (MPI). At every 
timestep the blocks that share border with other processors communicate their boundari es (ghost cel ls). Domain 
decomposition is made following a Hilbert ordering in similar fashion to the Ramses code ( Teyssier"|2 002 ) . Each block 



is assigned a number according to its position in a Hilbert curve of the rt\ evs — 1 order. Following such curve, the 
number of leaf blocks is divided evenly among the available processors, and the blocks are distributed to be solved. The 
tree structure is copied entirely to every processor, but only the leaf blocks are distributed. This, however, is a time 
consuming operation, and the work-load is only balanced every ni oa d timesteps (user defined typically ni oa d ~ 10 — 100), 
yielding an overhead of about 10 percent, depending on how fast the mesh structure changes on a particular problem. 

A. 5. Testing the code 

Several standard tests have been made to the co de, we present here one that w e believe is very illustrative: the classic 



double Mach reflection test problem proposed by |Woodward fe Colella (1984). The test was inspired on laboratory 
studies and consists of a planar shock that meets a reflecting surface a t an angle of tt/3. The reflec ting surface is taken 
to be the x— axis, and the initial conditions (explained in all detail in Woodward fc Colella||1984 ) of the problem are: 



(n 7J v P ) - j ( L4 > °> °> !) for x > x *, (Al6) 

W,n-\fa 8 .25sina, -8.25 cos a, 116.5) otherwise. ^ i0 ^ 

which correspond to a Mach 10 shock in air (using an adiabatic equation of state with C v = 2.5, 7 = 1.4). The 
problem is solved on a Cartesian domain [0, 4] x [0, 1] with the following boundary conditions: outflow boundaries 
in x = 4; reflective boundaries in y = for x > 1/6; at x = 0, and ?/ = 0forx<l/6 inflow conditions are used; 
at y = 1 time dependent inflow conditions are applied, tracing the position of the shock front that moves according 
to x s (t) = 1/6 + 10 y sin a + ?/tana. The flow configuration produces a series of shocks, two of which are almost 
perpendicular and are separated by a contact discontinuity, and a small jet at the base of the reflecting surface. This 
jet is driven by a pressure gradient and its evolution is known to be ve ry sensitive to numerical diffusion. The results 
of this test can be compared with those of Woodward fc Colella| ([1984]) , and many other codes which have adopted it 



as a standard test for two dimensional hydrodynamics (e.g. IMignone et al. 2007 



Stone et al.||2008| ). 

In Figure [6j we show the results of the double Mach reflection test with the HLLC solver, for an integration time 
of t = 0.2, for two different resolutions. Both runs start with 4 square root blocks aligned along the x-axis. The root 
blocks have 10 x 10 cells, covering a physical size of (1 x 1). The run in panel (a) has 5 levels of refinement, for a 
maximum resolution of Ax = Ay = 1/160; the run in panel (b) has 6 levels of refinement, and a maximum resolution 
of Ax = Ay = 1/320. Below each of the panels, the AMR structure (only the [0,3] x [0,1] region) is shown, each 
square corresponds to one n x x n y block, and in color we show the domain decomposition in 8 processors. Notice how 
the mesh is only refined at places where the flow properties change, and it remains coarse at regions where the flow is 
smooth. 

Another feature that can be identified is the fact that the jet formed along the reflecting surface (at x ~ 2.4-2.8) 
reaches the shock that is almost perpendicular to the surface (x ~ 2.8), bending it slightly at its base. This is often 
called "the kinked Mach stem" , and it is u nphysical, ca used by insufficient numerical diffusion at places where shocks 
closely align with the computational grid QQuirk|[l994| ). This issue becomes particularly critical for high resolution 
simulations (less numerical/artificial viscosity). 

A. 5.1. Parallel performance 

Is is customary to measure the performance of a parallel code in terms of its scaling with the number of processors. 
Ideally the computing time would scale as oc l/iV p , where N p is the number of processors. In fact, many fixed grid 
codes show near ideal scaling. The reason for such a good scaling is twofold: (i) most of the code is parallelizable, 
which restricts the communication between processors to passing boundaries and determining the timestep, and (ii) 
since the grid is fixed, each processor knows beforehand with which processors the communication will take place. 

Adaptive mesh codes use a significant amount of CPU time to maintain the AMR structure, and neighboring blocks 
can be in a different processors at different times. This causes AMR codes to scale less than ideally, but the CPU time 
and memory requirements can remain significantly lower than for fixed grid codes. Also, there are many variables that 
can impact the performance, and many of them are problem dependent. To illustrate this, we took the double Mach 
reflection test with an equivalent resolution of 2560 x 640 cells and ran it on 2, 4, 8, 16, and 32 processors. To achieve 
the same resolution at the highest level of refinement we used three combinations: 5 levels of refinement with 40 x 40 
cells blocks, 6 levels with blocks of 20 x 20 cells, and 7 levels with 10 x 10 blocks, the scaling results are presented in 
Figure [7| 
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Fig. 6. — Density contours and computational mesh at t = 0.2 for the double Mach reflection test, with the HLLC solver, at two different 
resolutions. The computational domain is [0,4] x [0,1], but for visual purposes we only display the [0,3] x [0,1] region. The equivalent 
resolution at the finest level is indicated in each panel. Below the contours we show the AMR blocks (each of 10 x 10 cells), and with colors 
we indicate the different processors in which the work load was distributed at the time. 
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Fig. 7. — Scaling of the Walicxe code for the double Mach reflection test at an equivalent resolution of 2560 x 640 cells, on 2, 4, 8, 16, 
and 32 processors. The different symbols are for a different combination of levels of refinement and block size, as indicated in the legend. 
We have included (dashed line) the ideal oc 1/N P scaling for reference. 



We see from the figure a good performance of the code, almost linear for a small number of processors. As we 
increase the number of processors the portion of domain that is solved by each processors becomes too small and the 
communication time becomes important compared to the time to solve each region. Certainly, if we require an even 
finer resolution, the transition from ideal, to not-so-ideal would happen at a larger number of processors. 

Is is quite noticeable that by allowing the mesh to be finer the CPU time drops dramatically. We must note that the 
outcome is particularly good because the problem is well suited for an adaptive mesh, the regions of sharp changes are 
well localized and are a small portion of the entire computational domain. If refinement is required in a larger portion 
of the domain, larger blocks and a small number of levels would be beneficial. 



A. 6. Fixing the kinked Mach stem ("carbuncle") problem 



The so-called numerical shock instability, or "carbuncle" probl em was first reported 



by|Quirk|(|1994[), a nd although 
& kitamura|2008p , researchers 
el" 



a full understanding of the instability and its cure is not reached ( Liou|2QQQ [Nishikawa < 

seem to agree that the culprit is insufficient viscosity when shocks align with the computational cell interfaces . Also, 
Riemann solvers that perform exceptionally well on shear layer problems (such as the Roe solver, |Roe||198l| ) suffer 
more of this pathology, while more diffusive ones, such as the HLL solver have been considered virtually carbuncle- free. 
In recent years, several remedies have been p ut forward to solve this instability, most of them for the solver of Roe 
(see for instance Nishikawa fc Kitamura|2QQ8 ). Among the remedies, a family of methods is to produce hybrid fluxes, 
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Fig. 8. — Density maps for the [2.25,2.75] x [0,0.5] region of the double Mach reflection test with the (a) HLLC, (b) the HLL, and (c) 
the hybrid HLL-HLLC Riemann solvers. The density values are indicated by the bar at the right. All the runs were done with a maximum 
resolution of Ax = Ay = 1/(640). 

combining a solver that is more diffusive, su ch as the HLL, with a more accurate, yet carbuncle-prone solver. In 



particular ( |Kim et al.|2009 Huang et al. 



2010), have proposed a combination of the HLL and HLLC fluxes. We have 



(A17) 



implemented a similar hybrid flux formed by a linear combination of the two fluxes 

Fhll-hllc = /^iFrllc + /^Frll, 

where the /3i + fa = 1, these coefficients are chosen in order to give more weight to the HLLC solver when shear flows 
are present, and to the HLL in the presence of strong shocks aligned with the grid. To do this we check the direction 
of the velocity difference vector between the left and right states: Aq = (ur — v>l,vr — vl). I n particular we use the 
magnitude of the projection of the unit vector in the direction of Aq with the unit vector normal to the cell interface, 
that is : oi\ = \ur — ul\/[{ur — ul) 2 + (vr — vl) 2 ] 1 ^ 2 - After some tests, we choose f3\ = y/a 1 , in figure pi we compare 
of the results from the double Mach reflection test with the HLLC, HLL, and HLL-HLLC solvers, all at t = 0.2, and 
all with the same equivalent resolution Ax = Ay = 1/ (640) . 

We can see from the Figure (panel a) how the HLLC solver alone produces a kinked Mach stem, but resolves sharply 
the contact discontinuity (the density change at an angle of 40deg. In fact some Kelvin-Helmholtz "rolls" start to 
appear, along this discontinuity and extend to the jet propagating along the reflecting surface (a:- axis). For the HLL 
solver (panel b) the kinked Mach stem is controlled (barely noticeable), but the diffusion at the contact surfaces is 
evident, producing a smeared solution compared to the more accurate HLLC solver. In panel (c) we show the results 
with the hybrid fluxes, and the best of the other two solvers is evident, a barely noticeable kinked Mach stem, with a 
sharp contact surface. 
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